All plotting from BLTrackfile¶
The file (from Fluka) has the format of
!BLTrackfile created by Ting
!x y z Px Py Pz t PDGid EvNum TrkId Parent weight
!cm cm cm GeV/c GeV/c GeV/c ns - - - - -
In [1]:
namelist = ['x', 'y', 'z', 'Px', 'Py', 'Pz', 't', 'PDGid', 'EventID', 'TrackID', 'ParentID', 'Weight']
units = ['cm', 'cm', 'cm', 'GeV/c', 'GeV/c', 'GeV/c', 'ns', '-', '-', '-', '-', '-']
In [2]:
# libraries
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as pc
from scipy.constants import c, R, k
from decimal import Decimal, getcontext
getcontext().prec = 16
import time
from scipy.stats import norm
In [3]:
# plotting help
from IPython.display import HTML
import io
import base64
class FlowLayout(object):
''' A class / object to display plots in a horizontal / flow layout below a cell '''
def __init__(self):
# string buffer for the HTML: initially some CSS; images to be appended
self.sHtml = """
<style>
.floating-box {
display: inline-block;
margin: 10px;
border: 3px solid #888888;
}
</style>
"""
def add_plot(self, oAxes):
''' Saves a PNG representation of a Matplotlib Axes object '''
Bio=io.BytesIO() # bytes buffer for the plot
fig = oAxes.get_figure()
fig.canvas.print_png(Bio) # make a png of the plot in the buffer
# encode the bytes as string using base 64
sB64Img = base64.b64encode(Bio.getvalue()).decode()
self.sHtml+= (
'<div class="floating-box">'+
'<img src="data:image/png;base64,{}\n">'.format(sB64Img)+
'</div>')
def PassHtmlToCell(self):
''' Final step - display the accumulated HTML '''
display(HTML(self.sHtml))
Load files¶
In [4]:
# number_of_files = np.array([1,7]).astype(int)
# number_of_files = ['+','-']
# number_of_files = ['Plus_1M']
number_of_files = [
'piplus_in_dt_B'
,
'piplus_in_dt_noB'
]
total_number_of_files = len(number_of_files)
In [5]:
data = []
for number in number_of_files:
# # number = int(number)
# filename = f'{name_of_files}_{str(number)}.txt' # input your file name
# for name in name_of_files:
filename = f'{number}.txt' # input your file name
print(f' loading "{filename}" ...')
t0 = time.time()
data.append(np.loadtxt(filename))
print('This file has' , np.loadtxt(filename).shape[0], 'number of pions')
t1 = time.time()
t = int(t1-t0)
print(f' loaded "{filename}" ..., it took:{t}s')
# cut to minimum row
min_line = []
for d in data:
min_line.append(len(d))
line = min(min_line)
data_np = []
for i in range(total_number_of_files):
data_np.append( data[i][0:line,:] )
data = np.array(data_np)
del data_np
data.shape
loading "piplus_in_dt_B.txt" ... This file has 1072572 number of pions loaded "piplus_in_dt_B.txt" ..., it took:15s loading "piplus_in_dt_noB.txt" ... This file has 799707 number of pions loaded "piplus_in_dt_noB.txt" ..., it took:11s
Out[5]:
(2, 799707, 12)
Define variable and collect data needed¶
In [6]:
all_data = np.empty([total_number_of_files, line, 9])
x, y, r, emission_angle, Px, Py, Pz, P, PDGid = all_data[:,:,0], all_data[:,:,1], all_data[:,:,2], all_data[:,:,3], all_data[:,:,4], all_data[:,:,5], all_data[:,:,6], all_data[:,:,7], all_data[:,:,8]
for i in range(total_number_of_files):
x[i], y[i] = data[i,:,0], data[i,:,1]
r[i] = np.sqrt(x[i]**2 + y[i]**2)
Px[i], Py[i], Pz[i]= data[i,:,3], data[i,:,4], data[i,:,5]
P[i] = np.sqrt(Pz[i]**2+Px[i]**2+Py[i]**2)
emission_angle[i] = np.arctan( np.sqrt(Px[i]**2+Py[i]**2)/Pz[i] )
PDGid[i] = data[i,:,7]
In [7]:
# parameters same for all data
nu_e = PDGid == 12
mu_plus = PDGid == -13
nu_mu_bar = PDGid == -14
pion_plus = PDGid == 211
pion_minus = PDGid == -211
In [8]:
Px.shape
Out[8]:
(2, 799707)
Define a cut here¶
In [9]:
# multiple cut
cut=[]
pion720880 = ((P>0.720)&(P<0.880)&pion_plus)
emit_xy = ( (x>-10)&(x<10)&(Px/Pz<15/1000)&(Px/Pz>-15/1000)& (y>-10)&(y<10)&(Py/Pz<15/1000)&(Py/Pz>-15/1000) )
cut_list = [pion_plus, pion720880, emit_xy] # [ nu_e , nu_mu_bar, nu_e006, nu_mu_bar006] , emit_xy
cut_name_list = [r'$\pi^+$ ',r'$\pi^+$ 720<P<880', 'emit_xy'] # , 'emit_xy'
total_number_of_cut = len(cut_name_list)
for c in cut_list:
cut.append(c)
Selection of data of which file¶
In [10]:
x_s_c, y_s_c, r_s_c, emission_angle_s_c, Px_s_c, Py_s_c, Pz_s_c, P_s_c, PDGid_s_c = [], [], [], [], [], [], [], [], []
for cut in cut_list:
selected_data=[]
for i in range(total_number_of_files):
selected_data.append( all_data[i][cut[i]] )
x_s, y_s, r_s, emission_angle_s, Px_s, Py_s, Pz_s, P_s, PDGid_s = [], [], [], [], [], [], [], [], []
for i in range(total_number_of_files):
x_s.append(selected_data[i][:,0])
y_s.append(selected_data[i][:,1])
r_s.append(selected_data[i][:,2])
emission_angle_s.append(selected_data[i][:,3])
Px_s.append(selected_data[i][:,4])
Py_s.append(selected_data[i][:,5])
Pz_s.append(selected_data[i][:,6])
P_s.append(selected_data[i][:,7])
PDGid_s.append(selected_data[i][:,8])
x_s_c.append(x_s)
y_s_c.append(y_s)
r_s_c.append(r_s)
emission_angle_s_c.append(emission_angle_s)
Px_s_c.append(Px_s)
Py_s_c.append(Py_s)
Pz_s_c.append(Pz_s)
P_s_c.append(P_s)
PDGid_s_c.append(PDGid_s)
In [11]:
namelist_s_c = ['x', 'y', 'r', r'$\theta$','Px', 'Py', 'Pz', 'P', 'PDGid']
len(namelist_s_c)
Out[11]:
9
In [12]:
All_data_s_c = [x_s_c, y_s_c, r_s_c, emission_angle_s_c, Px_s_c, Py_s_c, Pz_s_c, P_s_c, PDGid_s_c ]
Plotting¶
Momentum distribution¶
In [13]:
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k','pink','brown','olive','orange', 'purple']
linestyle = "-"
bin_min = 0
bin_max = 2000
bin_width = 10
bins_number = int(bin_max/bin_width+1)
bins_number
bin_P = np.linspace(bin_min, bin_max, bins_number)
In [14]:
total_number_of_files
Out[14]:
2
In [15]:
fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(111)
# bin = 100
factor = 1000
## plot type
for i in range(total_number_of_files):
for ic in range(total_number_of_cut):
# if ic == 2:
# linestyle = ':'
# n, bins, patches = ax.hist(P_s_c[ic][i], bins=bin, histtype='step', color = colors[ic-int(len(cut_list)/2)], linestyle=linestyle,
# label=fr'{cut_name_list[ic]} : Total number is {P_s_c[ic][i].shape[0]}')
# else:
# linestyle = '-'
n, bins, patches = ax.hist(P_s_c[ic][i]*factor, bins=bin_P, histtype='step', color = colors[total_number_of_cut*i +ic], linestyle=linestyle,
label=fr'{number_of_files[i]}_{cut_name_list[ic]} : Total number is {P_s_c[ic][i].shape[0]}')
## optional config
plt.legend(loc='upper right')
plt.xlabel('P(MeV)')
plt.ylabel(f' N/{bin_width} MeV')
# ax.set_ylim(0,1000)
# plt.title(f'{cut_name_list[i]}') #f'$10^{int(np.log10(data.shape[1]))}$ {number_of_files[i]}00MeV muon decays. ~14% in the narrow cone.')
save_name = f'{data.shape[1]}_pion_momentum_dis.png'
plt.savefig(save_name)
Emission angle distribution¶
In [16]:
bin_min = 0
bin_max = np.pi/2
bin_width = 0.02
bins_number = int(bin_max/bin_width+1)
bins_number
bin_theta = np.linspace(bin_min, bin_max, bins_number)
print(bin_theta.shape)
bin_theta
(79,)
Out[16]:
array([0. , 0.02013841, 0.04027683, 0.06041524, 0.08055366,
0.10069207, 0.12083049, 0.1409689 , 0.16110732, 0.18124573,
0.20138414, 0.22152256, 0.24166097, 0.26179939, 0.2819378 ,
0.30207622, 0.32221463, 0.34235305, 0.36249146, 0.38262987,
0.40276829, 0.4229067 , 0.44304512, 0.46318353, 0.48332195,
0.50346036, 0.52359878, 0.54373719, 0.5638756 , 0.58401402,
0.60415243, 0.62429085, 0.64442926, 0.66456768, 0.68470609,
0.70484451, 0.72498292, 0.74512133, 0.76525975, 0.78539816,
0.80553658, 0.82567499, 0.84581341, 0.86595182, 0.88609024,
0.90622865, 0.92636706, 0.94650548, 0.96664389, 0.98678231,
1.00692072, 1.02705914, 1.04719755, 1.06733597, 1.08747438,
1.10761279, 1.12775121, 1.14788962, 1.16802804, 1.18816645,
1.20830487, 1.22844328, 1.2485817 , 1.26872011, 1.28885852,
1.30899694, 1.32913535, 1.34927377, 1.36941218, 1.3895506 ,
1.40968901, 1.42982743, 1.44996584, 1.47010425, 1.49024267,
1.51038108, 1.5305195 , 1.55065791, 1.57079633])
In [17]:
fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(111)
## plot type
for i in range(total_number_of_files):
for ic in range(len(cut_list)):
n, bins, patches = ax.hist(emission_angle_s_c[ic][i], bins=bin_theta, histtype='step', color = colors[i*total_number_of_cut + ic], linestyle=linestyle,
label=fr'{number_of_files[i]}_{cut_name_list[ic]} : Total number is {emission_angle_s_c[ic][i].shape[0]}')
# if ic ==1:
# print(r'number of pion in $\theta<0.1$ is', sum(n[bins[1::]<0.1]))
## optional config
ax.legend()
# plt.legend(legend_labels, loc='upper right')
plt.xlabel(fr'$\theta_{{lab}}$ (rads)')
plt.ylabel(fr' N / {bin_width} rad')
plt.title(f'{cut_name_list[i]}')
save_name = f'{data.shape[1]}_emission_angle_dis.png'
# plt.yscale("log")
# save_name = f'Emission_angle_distribution_{data.shape[1]} decay.png'
plt.savefig(save_name)
2D plot¶
P-theta Heat map¶
In [18]:
# functioni call
def get_histogram2d(x=None, y=None, z=None,
bins=10, range=None,
xscale="linear", yscale="linear", cscale="linear",
normed=False, cmap=None, clim=(None, None),
ax1=None, grid=True, shading='flat', colorbar={},
cbi_kwargs={'orientation': 'vertical'},
xlabel="", ylabel="", clabel="", title="",
fname="hist2d.png"):
"""
creates a 2d histogram
Parameters
----------
x, y, z :
x and y coordinaten for z value, if z is None the 2d histogram of x and z is calculated
numpy.histogram2d parameters:
range : array_like, shape(2,2), optional
bins : int or array_like or [int, int] or [array, array], optional
ax1: mplt.axes
if None (default) a olt.figure is created and histogram is stored
if axis is give, the axis and a pcolormesh object is returned
colorbar : dict
plt.pcolormesh parameters:
clim=(vmin, vmax) : scalar, optional, default: clim=(None, None)
shading : {'flat', 'gouraud'}, optional
normed: string
colum, row, colum1, row1 (default: None)
{x,y,c}scale: string
'linear', 'log' (default: 'linear')
"""
if z is None and (x is None or y is None):
sys.exit("z and (x or y) are all None")
if ax1 is None:
fig, ax = plt.subplots(1)
fig.subplots_adjust(wspace=0.3, left=0.05, right=0.95)
else:
ax = ax1
if z is None:
z, xedges, yedges = np.histogram2d(x, y, bins=bins, range=range)
z = z.T
else:
xedges, yedges = x, y
if normed:
if normed == "colum":
z = z / np.sum(z, axis=0)
elif normed == "row":
z = z / np.sum(z, axis=1)[:, None]
elif normed == "colum1":
z = z / np.amax(z, axis=0)
elif normed == "row1":
z = z / np.amax(z, axis=1)[:, None]
else:
sys.exit("Normalisation %s is not known.")
color_norm = mpl.colors.LogNorm() if cscale == "log" else None
vmin, vmax = clim
im = ax.pcolormesh(xedges, yedges, z, shading=shading, vmin=vmin, vmax=vmax, norm=color_norm, cmap=cmap)
if colorbar is not None:
cbi = plt.colorbar(im, **cbi_kwargs)
cbi.ax.tick_params(axis='both', **{"labelsize": 14})
cbi.set_label(clabel)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_xscale(xscale)
ax.set_yscale(yscale)
ax.set_title(title)
return fig, ax, im
In [19]:
bin_min = 0
# bin_max = max(number_of_files)*100
bin_max = 2000
bin_width = 10
bins_number = int(bin_max/bin_width+1)
bin_P = np.linspace(bin_min, bin_max, bins_number)
bin_min = 0
bin_max = np.pi/2
bin_width = 0.02
bins_number = int(bin_max/bin_width+1)
bin_theta = np.linspace(bin_min, bin_max, bins_number)
bins_P_theta = [bin_P, bin_theta]
In [20]:
# main plotting
oPlot = FlowLayout() # create an empty FlowLayout
#Heat map for r theta
#####################
for i in range(total_number_of_files):
ic = 0
run_name = fr'{i}'
plot_title = fr"Heatmap of P-$\theta$ of {number_of_files[i]}"
xlabel = f"P (MeV)"
ylabel = fr"$\theta$ (rad)"
cmap = "Blues"# "coolwarm" # "BuPu" "YlOrRd"
# bins = 100
cscale = "linear"
file_name = f"plots/scatter_2dhistogram_{run_name}_cscale{cscale}.png"
# Also plot a heatmap of the scatter plot instead of just dots
fig, ax, im = get_histogram2d(P_s_c[ic][i]*factor, emission_angle_s_c[ic][i],
title=plot_title, xlabel=xlabel, ylabel=ylabel, bins=bins_P_theta,
cmap=cmap, cscale=cscale)
# # Plot a black line through the middle
# xmin = min(shower_energy_log10)
# xmax = max(shower_energy_log10)
# ymin = min(shower_energy_log10_predict)
# ymax = max(shower_energy_log10_predict)
# ax.plot([min(xmin, ymin), max(xmax, ymax)], [min(xmin, ymin), max(xmax, ymax)], 'k--')
# ax.set_ylim(0,np.pi/2)
plt.xlim(0,1000)
plt.ylim(0,1)
plt.tight_layout()
save_name = f'{data.shape[1]}_{number_of_files[i]}_P_theta_2Ddis.png'
plt.savefig(save_name, dpi=300)
oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
plt.close() # this gets rid of the plot so it doesn't appear in the cell
oPlot.PassHtmlToCell()
In [21]:
oPlot = FlowLayout() # create an empty FlowLayout
# fig = plt.figure(figsize=(10,10))
# ax = fig.add_subplot(121)
## plot type
for i in range(total_number_of_files):
for k in range(3):
fig, ax = plt.subplots(1, 1, figsize=(6,6))
for ic in range(total_number_of_cut):
if k ==2:
ax.scatter(r_s_c[ic][i], emission_angle_s_c[ic][i], color = colors[ic], s=0.1, label=fr'{cut_name_list[ic]} : Total number is {P_s_c[ic][i].shape[0]}')
ax.set_ylabel(f'{namelist_s_c[k+1]}(rad)')
else:
ax.scatter(All_data_s_c[k][ic][i], All_data_s_c[k+4][ic][i]/Pz_s_c[ic][i], color = colors[ic], s=0.1, label=fr'{cut_name_list[ic]} : Total number is {P_s_c[ic][i].shape[0]}')
ax.set_ylabel(f'{namelist_s_c[k]}\'(rad)')
## optional config
ax.legend(loc='upper right')
ax.set_xlabel(f'{namelist_s_c[k]} (cm)')
plt.title(f'{number_of_files[i]}')
save_name = f'{data.shape[1]}_{number_of_files[i]}_emittance_{namelist_s_c[k]}.png'
plt.savefig(save_name)
oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
plt.close() # this gets rid of the plot so it doesn't appear in the cell
oPlot.PassHtmlToCell()
In [22]:
# oPlot = FlowLayout() # create an empty FlowLayout
# i = 0
# # fig = plt.figure(figsize=(10,10))
# # ax = fig.add_subplot(121)
# ## plot type
# fig, ax = plt.subplots(1, 1, figsize=(7,7))
# for ic in range(len(cut_list)):
# ax.scatter(r_s_c[ic][i], emission_angle_s_c[ic][i], color = colors[ic], s=0.1, label=fr'{cut_name_list[ic]} : Total number is {P_s_c[ic][i].shape[0]}')
# ## optional config
# ax.legend(loc='upper right')
# ax.set_xlabel(f'r (cm)')
# ax.set_ylabel(r'$\theta$ (rad)')
# save_name = f'{data.shape[1]} {name_of_files}{number_of_files[i]}_emittance_r_{namelist_s_c[k]}.png'
# plt.savefig(save_name)
# oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
# plt.close() # this gets rid of the plot so it doesn't appear in the cell
# oPlot.PassHtmlToCell()
In [23]:
bin_min = -200
# bin_max = max(number_of_files)*100
bin_max = 200
bin_width_xy = 1
bins_number = int(bin_max/bin_width_xy+1)
bin_xy = np.linspace(bin_min, bin_max, bins_number)
bin_min = -1000
bin_max = 1000
bin_width_theta = 1
bins_number = int(bin_max/bin_width_theta+1)
bin_theta = np.linspace(bin_min, bin_max, bins_number)
bins_xy_theta = [bin_xy, bin_theta]
Emittance and confidence level fit¶
https://matplotlib.org/stable/gallery/statistics/confidence_ellipse.html
In [24]:
from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms
def confidence_ellipse(x, y, ax, n_std=3.0, facecolor='none', **kwargs):
"""
Create a plot of the covariance confidence ellipse of `x` and `y`
Parameters
----------
x, y : array_like, shape (n, )
Input data.
ax : matplotlib.axes.Axes
The axes object to draw the ellipse into.
n_std : float
The number of standard deviations to determine the ellipse's radiuses.
Returns
-------
matplotlib.patches.Ellipse
Other parameters
----------------
kwargs : `~matplotlib.patches.Patch` properties
"""
if x.size != y.size:
raise ValueError("x and y must be the same size")
cov = np.cov(x, y)
pearson = cov[0, 1]/np.sqrt(cov[0, 0] * cov[1, 1])
# Using a special case to obtain the eigenvalues of this
# two-dimensionl dataset.
ell_radius_x = np.sqrt(1 + pearson)
ell_radius_y = np.sqrt(1 - pearson)
ellipse = Ellipse((0, 0),
width=ell_radius_x * 2,
height=ell_radius_y * 2,
facecolor=facecolor,
**kwargs)
# Calculating the stdandard deviation of x from
# the squareroot of the variance and multiplying
# with the given number of standard deviations.
scale_x = np.sqrt(cov[0, 0]) * n_std
mean_x = np.mean(x)
# calculating the stdandard deviation of y ...
scale_y = np.sqrt(cov[1, 1]) * n_std
mean_y = np.mean(y)
transf = transforms.Affine2D() \
.rotate_deg(45) \
.scale(scale_x, scale_y) \
.translate(mean_x, mean_y)
ellipse.set_transform(transf + ax.transData)
return ax.add_patch(ellipse)
In [25]:
# main plotting
oPlot = FlowLayout() # create an empty FlowLayout
#Heat map for r theta
#####################
print(" All momentum cut emittacnce cut")
for i in range(total_number_of_files):
print(number_of_files[i])
for k in range(3):
print(" ",namelist_s_c[k])
for ic in range(len(cut_list)):
run_name = fr' {cut_name_list[ic]}'
plot_title = fr"Heatmap of emittance of {run_name}, bin width = {bin_width_xy}cm*{bin_width_theta}mrad." +"\n" +f"Total number is {x_s_c[ic][i].shape[0]}"
xlabel = f"{namelist_s_c[k]} (cm)"
ylabel = f"{namelist_s_c[k]}\'(mrad)"
cmap = "Reds"#"coolwarm" # "BuPu" "YlOrRd"
# bins = 100
cscale = "linear"
file_name = f"plots/scatter_2dhistogram_{run_name}_cscale{cscale}.png"
# Also plot a heatmap of the scatter plot instead of just dots
if k ==2:
x = r_s_c[ic][i]
y = emission_angle_s_c[ic][i]*1e3
fig, ax, im = get_histogram2d(x, y, bins=bins_xy_theta, title=plot_title, xlabel=xlabel, ylabel=ylabel, cmap=cmap, cscale=cscale)
confidence_ellipse(x, y, ax, n_std=1, label=r'$1\sigma$', edgecolor='green')
ax.set_ylabel(f'{namelist_s_c[k+1]}(mrad)')
else:
x = All_data_s_c[k][ic][i]
y = All_data_s_c[k+4][ic][i]/Pz_s_c[ic][i]*1e3
fig, ax, im = get_histogram2d(x, y, bins=bins_xy_theta, title=plot_title, xlabel=xlabel, ylabel=ylabel, cmap=cmap, cscale=cscale)
confidence_ellipse(x, y, ax, n_std=1, label=r'$1\sigma$', edgecolor='green')
ax.set_ylabel(f'{namelist_s_c[k]}\'(mrad)')
# # Plot a black line through the middle
# xmin = min(shower_energy_log10)
# xmax = max(shower_energy_log10)
# ymin = min(shower_energy_log10_predict)
# ymax = max(shower_energy_log10_predict)
# ax.plot([min(xmin, ymin), max(xmax, ymax)], [min(xmin, ymin), max(xmax, ymax)], 'k--')
# ax.set_ylim(0,np.pi/2)
# plt.xlim(-50, 50)
# plt.ylim(-50, 50)
plt.legend()
plt.tight_layout()
save_name = f'{data.shape[1]} {number_of_files[i]}_{namelist_s_c[k]}_{plot_title}.png'
# plt.savefig(save_name, dpi=300)
oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
plt.close() # this gets rid of the plot so it doesn't appear in the cell
oPlot.PassHtmlToCell()
All momentum cut emittacnce cut
piplus_in_dt_B
x
y
r
piplus_in_dt_noB
x
y
r
Emittance and Acceptance fit¶
In [26]:
# main plotting
oPlot = FlowLayout() # create an empty FlowLayout
#Heat map for r theta
#####################
# fit parameter
u=0 #x-position of the center
v=0 #y-position of the center
a=10 #radius on the x-axis
b=15 #radius on the y-axis
t_rot=-15*np.pi/180 #rotation angle
t = np.linspace(0, 2*np.pi, 100)
Ell = np.array([a*np.cos(t) , b*np.sin(t)])
#u,v removed to keep the same center location
R_rot = np.array([[np.cos(t_rot) , -np.sin(t_rot)],[np.sin(t_rot) , np.cos(t_rot)]])
#2-D rotation matrix
Ell_rot = np.zeros((2,Ell.shape[1]))
for k in range(Ell.shape[1]):
Ell_rot[:,k] = np.dot(R_rot,Ell[:,k])
# main plot
print(" All momentum cut emittacnce cut")
for i in range(total_number_of_files):
print(number_of_files[i])
for k in range(3):
print(" ",namelist_s_c[k])
for ic in range(len(cut_list)):
run_name = fr' {cut_name_list[ic]}'
plot_title = fr"Heatmap of emittance of {run_name}, bin width = {bin_width_xy}cm*{bin_width_theta}mrad." +"\n" +f"Total number is {x_s_c[ic][i].shape[0]}"
xlabel = f"{namelist_s_c[k]} (cm)"
ylabel = f"{namelist_s_c[k]}\'(mrad)"
cmap = "Reds"#"coolwarm" # "BuPu" "YlOrRd"
# bins = 100
cscale = "linear"
file_name = f"plots/scatter_2dhistogram_{run_name}_cscale{cscale}.png"
# Also plot a heatmap of the scatter plot instead of just dots
if k ==2:
x = r_s_c[ic][i]
y = emission_angle_s_c[ic][i]*1e3
fig, ax, im = get_histogram2d(x, y, bins=bins_xy_theta, title=plot_title, xlabel=xlabel, ylabel=ylabel, cmap=cmap, cscale=cscale)
ax.set_ylabel(f'{namelist_s_c[k+1]}(mrad)')
plt.xlim(0, 100)
plt.ylim(0, 250)
else:
x = All_data_s_c[k][ic][i]
y = All_data_s_c[k+4][ic][i]/Pz_s_c[ic][i]*1e3
fig, ax, im = get_histogram2d(x, y, bins=bins_xy_theta, title=plot_title, xlabel=xlabel, ylabel=ylabel, cmap=cmap, cscale=cscale)
ax.set_ylabel(f'{namelist_s_c[k]}\'(mrad)')
plt.plot( u+Ell_rot[0,:] , v+Ell_rot[1,:],'blue', label=fr'area = {a}cm*{b}mrad*$\pi$ ~ 1.5 mmrad' ) #rotated ellipse darkorange
plt.xlim(-50, 50)
plt.ylim(-50, 50)
# plt.plot( u+Ell[0,:] , v+Ell[1,:], color = 'g' ) #initial ellipse
plt.legend()
plt.tight_layout()
save_name = f'{data.shape[1]} {number_of_files[i]}_{namelist_s_c[k]}_{plot_title}.png'
# plt.savefig(save_name, dpi=300)
oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
plt.close() # this gets rid of the plot so it doesn't appear in the cell
oPlot.PassHtmlToCell()
All momentum cut emittacnce cut
piplus_in_dt_B
x
y
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
r
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
piplus_in_dt_noB
x
y
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
r
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
heat map xy¶
In [27]:
bin_min = -200
# bin_max = max(number_of_files)*100
bin_max = 200
bin_width_xy = 1
bins_number = int(bin_max/bin_width_xy+1)
bin_xy = np.linspace(bin_min, bin_max, bins_number)
bins_xy_xy = [bin_xy, bin_xy]
# bins = 50
In [28]:
# main plotting
oPlot = FlowLayout() # create an empty FlowLayout
#Heat map for r theta
#####################
print(" All momentum cut emittacnce cut")
for i in range(total_number_of_files):
print(number_of_files[i])
for k in range(2):
print(" ",namelist_s_c[k])
for ic in range(len(cut_list)):
run_name = fr' {cut_name_list[ic]}'
plot_title = fr"Heatmap of emittance of {run_name}, bin width = {bin_width_xy}cm*{bin_width_theta}mrad." +"\n" +f"Total number is {x_s_c[ic][i].shape[0]}"
xlabel = f"x (cm)"
ylabel = f"y'(cm)"
cmap = "Greens"#"coolwarm" # "BuPu" "YlOrRd"
cscale = "linear"
file_name = f"plots/scatter_2dhistogram_{run_name}_cscale{cscale}.png"
# Also plot a heatmap of the scatter plot instead of just dots
x, y = x_s_c[ic][i], y_s_c[ic][i]
fig, ax, im = get_histogram2d(x, y, bins=bins_xy_xy, title=plot_title, xlabel=xlabel, ylabel=ylabel, cmap=cmap, cscale=cscale)
confidence_ellipse(x, y, ax, n_std=1, label=r'$1\sigma$', edgecolor='red')
# xy_range = 20
# plt.xlim(-xy_range, xy_range)
# plt.ylim(-xy_range, xy_range)
plt.tight_layout()
save_name = f'{data.shape[1]} {number_of_files[i]}_xy_{plot_title}.png'
# plt.savefig(save_name, dpi=300)
oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
plt.close() # this gets rid of the plot so it doesn't appear in the cell
oPlot.PassHtmlToCell()
All momentum cut emittacnce cut
piplus_in_dt_B
x
y
piplus_in_dt_noB
x
y
In [29]:
# main plotting
oPlot = FlowLayout() # create an empty FlowLayout
#Heat map for r theta
#####################
print(" All momentum cut emittacnce cut")
for i in range(total_number_of_files):
print(number_of_files[i])
for k in range(2):
print(" ",namelist_s_c[k])
for ic in range(len(cut_list)):
run_name = fr' {cut_name_list[ic]}'
plot_title = fr"Heatmap of emittance of {run_name}, bin width = {bin_width_xy}cm*{bin_width_theta}mrad." +"\n" +f"Total number is {x_s_c[ic][i].shape[0]}"
xlabel = f"x (cm)"
ylabel = f"y'(cm)"
cmap = "Greens"#"coolwarm" # "BuPu" "YlOrRd"
cscale = "linear"
file_name = f"plots/scatter_2dhistogram_{run_name}_cscale{cscale}.png"
# Also plot a heatmap of the scatter plot instead of just dots
x, y = x_s_c[ic][i], y_s_c[ic][i]
fig, ax, im = get_histogram2d(x, y, bins=bins_xy_xy, title=plot_title, xlabel=xlabel, ylabel=ylabel, cmap=cmap, cscale=cscale)
confidence_ellipse(x, y, ax, n_std=1, label=r'$1\sigma$', edgecolor='red')
xy_range = 20
plt.xlim(-xy_range, xy_range)
plt.ylim(-xy_range, xy_range)
plt.legend()
plt.tight_layout()
# save_name = f'{data.shape[1]} {number_of_files[i]}_{namelist_s_c[k]}_{plot_title}.png'
# plt.savefig(save_name, dpi=300)
oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
plt.close() # this gets rid of the plot so it doesn't appear in the cell
oPlot.PassHtmlToCell()
All momentum cut emittacnce cut
piplus_in_dt_B
x
y
piplus_in_dt_noB
x
y
x distribution¶
In [30]:
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111)
# bin = 100
## plot type
for i in range(total_number_of_files):
for ic in range(total_number_of_cut):
# if ic == 2:
# linestyle = ':'
# n, bins, patches = ax.hist(P_s_c[ic][i], bins=bin, histtype='step', color = colors[ic-int(len(cut_list)/2)], linestyle=linestyle,
# label=fr'{cut_name_list[ic]} : Total number is {P_s_c[ic][i].shape[0]}')
# else:
# linestyle = '-'
n, bins, patches = ax.hist(x_s_c[ic][i], bins=bin_xy, histtype='step', color = colors[total_number_of_cut*i +ic], linestyle=linestyle,
label=fr'{number_of_files[i]}_{cut_name_list[ic]} : Total number is {x_s_c[ic][i].shape[0]}')
ci = norm(*norm.fit(x_s_c[ic][i])).interval(0.68) # fit a normal distribution and get 95% c.i.
plt.fill_betweenx([0, n.max()], ci[0], ci[1], color=colors[total_number_of_cut*i +ic], alpha=0.1) # Mark between 0 and the highest bar in the histogram
## optional config
plt.legend(loc='upper right')
plt.xlabel('x (cm)')
plt.ylabel(f' N/{bin_width} MeV')
# ax.set_ylim(0,1000)
# plt.title(f'{cut_name_list[i]}') #f'$10^{int(np.log10(data.shape[1]))}$ {number_of_files[i]}00MeV muon decays. ~14% in the narrow cone.')
save_name = f'{data.shape[1]}_x_dis.png'
plt.savefig(save_name)
3D plots¶
In [39]:
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
# %matplotlib
%matplotlib widget
In [32]:
# plt.close("all")
In [40]:
# main plotting
# oPlot = FlowLayout() # create an empty FlowLayout
#Heat map for r theta
#####################
print(" All momentum cut emittacnce cut")
for i in range(total_number_of_files):
print(number_of_files[i])
fig = plt.figure(figsize=(18,12))
for k in range(2):
print(" ",namelist_s_c[k])
for ic in range(len(cut_list)):
run_name = fr' {cut_name_list[ic]}'
plot_title = fr"Heatmap of emittance of {run_name}, bin width = {bin_width_xy}cm*{bin_width_theta}mrad." +"\n" +f"Total number is {x_s_c[ic][i].shape[0]}"
xlabel = f"x (cm)"
ylabel = f"y'(cm)"
cmap = "Greens"#"coolwarm" # "BuPu" "YlOrRd"
cscale = "linear"
file_name = f"plots/scatter_2dhistogram_{run_name}_cscale{cscale}.png"
ax = plt.subplot(2,3,1+3*k+ic, projection='3d')
# Also plot a heatmap of the scatter plot instead of just dots
x, y = x_s_c[ic][i], y_s_c[ic][i]
# Make data
hist, xedges, yedges = np.histogram2d(x, y, bins=50)
X, Y = np.meshgrid(xedges[:-1] + 0.25, yedges[:-1] + 0.25)
Z = hist
# print(X.shape, Y.shape, Z.shape)
# Plot the surface
plt.style.use('_mpl-gallery')
cmap = plt.cm.get_cmap('viridis') # You can choose any colormap here
# Creating the figure
# fig = plt.figure(figsize=(5,5))
# ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, cmap=cm.Blues, edgecolor='none', alpha=0.7) #
# save_name = f'{data.shape[1]} {number_of_files[i]}_{namelist_s_c[k]}_{plot_title}.png'
# plt.savefig(save_name, dpi=300)
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.title(plot_title)
# Adding colorbar
# fig.colorbar(surf, ax=ax, label='Frequency')
# plt.show()
# oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
# plt.close() # this gets rid of the plot so it doesn't appear in the cell
plt.tight_layout()
# oPlot.PassHtmlToCell()
All momentum cut emittacnce cut
piplus_in_dt_B
x
y
piplus_in_dt_noB
x
y
In [34]:
# # main plotting
# # plt.rcParams.update({'font.size': 5})
# # oPlot = FlowLayout() # create an empty FlowLayout
# #Heat map for r theta
# #####################
# %matplotlib
# print(" All momentum cut emittacnce cut")
# for i in range(total_number_of_files):
# print(number_of_files[i])
# for k in range(2):
# print(" ",namelist_s_c[k])
# for ic in range(len(cut_list)):
# run_name = fr' {cut_name_list[ic]}'
# plot_title = fr"Heatmap of emittance of {run_name}, bin width = {bin_width_xy}cm*{bin_width_theta}mrad." +"\n" +f"Total number is {x_s_c[ic][i].shape[0]}"
# xlabel = f"x (cm)"
# ylabel = f"y'(cm)"
# cmap = "Greens"#"coolwarm" # "BuPu" "YlOrRd"
# cscale = "linear"
# file_name = f"plots/scatter_2dhistogram_{run_name}_cscale{cscale}.png"
# fig = plt.figure(figsize=(3,3))
# cmap = plt.cm.get_cmap('viridis') # You can choose any colormap here
# plt.style.use('_mpl-gallery')
# ax = fig.add_subplot(111, projection='3d')
# # Also plot a heatmap of the scatter plot instead of just dots
# if k ==2:
# x = r_s_c[ic][i]
# y = emission_angle_s_c[ic][i]*1e3
# # Make data
# hist, xedges, yedges = np.histogram2d(x, y, bins=50)
# X, Y = np.meshgrid(xedges[:-1] + 0.25, yedges[:-1] + 0.25)
# Z = hist
# surf = ax.plot_surface(X, Y, Z, cmap=cm.Blues, edgecolor='none', alpha=0.7) #
# ax.set_ylabel(f'{namelist_s_c[k+1]}(mrad)')
# else:
# x = All_data_s_c[k][ic][i]
# y = All_data_s_c[k+4][ic][i]/Pz_s_c[ic][i]*1e3
# # Make data
# hist, xedges, yedges = np.histogram2d(x, y, bins=50)
# X, Y = np.meshgrid(xedges[:-1] + 0.25, yedges[:-1] + 0.25)
# Z = hist
# surf = ax.plot_surface(X, Y, Z, cmap=cm.Blues, edgecolor='none', alpha=0.7) #
# ax.set_ylabel(f'{namelist_s_c[k]}\'(mrad)')
# # save_name = f'{data.shape[1]} {number_of_files[i]}_{namelist_s_c[k]}_{plot_title}.png'
# # plt.savefig(save_name, dpi=300)
# plt.xlabel(xlabel)
# plt.title(plot_title)
# # Adding colorbar
# # fig.colorbar(surf, ax=ax, label='Frequency', shrink=0.5)
# plt.tight_layout()
# oPlot.add_plot(ax) # pass it to the FlowLayout to save as an image
# plt.close() # this gets rid of the plot so it doesn't appear in the cell
# # oPlot.PassHtmlToCell()
In [35]:
!jupyter nbconvert --to html Plot_BL_files
[NbConvertApp] Converting notebook Plot_BL_files.ipynb to html [NbConvertApp] Writing 7827148 bytes to Plot_BL_files.html